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Abstract 

We propose a novel method for fitting planar B-spline curves 
to unorganized data points. In traditional methods, optimiza- 
tion of control points and foot points are performed in two very 
time-consuming steps in each iteration: 1) control points are 
updated by setting up and solving a linear system of equations; 
and 2) foot points are computed by projecting each data point 
onto a B-spline curve. Our method uses the L-BFGS optimiza- 
tion method to optimize control points and foot points simulta- 
neously and therefore it does not need to perform either matrix 
computation or foot point projection in every iteration. As a 
result, our method is much faster than existing methods. 



1 Introduction 

Curve fitting is a fundamental problem in many fields, such 
as computer graphics, image processing, shape modeling 
and data mining. Depending on applications, different 
types of curves such as parametric curves, implicit curves 
and subdivision curves are used for fitting. In this paper, 
we study the problem of fitting planar B-spline curves to 
unorganized data points. 

Given a set of unorganized data points {Xijflj C 
R 2 sampled from the outline of a planar shape, the 
aim of curve fitting is to find a B-spline curve P(t) = 
X/iLi PiNi(t) that best approximates the shape's outline. 
The outline is called a target shape, and the B-spline curve 
is called a fitting curve. Here, V := {Pi}" =1 C R 2 is the 
set of B-spline control points, {Ni(t)}" =1 are B-spline basis 
functions. We suppose that knots of the B-spline curve are 
fixed and therefore not subject to optimization, and all the 
basis functions are thus defined on fixed, uniform spaced 
knots throughout the curve fitting process. 

For a data point Xfc, let P(^) denote the nearest point 
of Xfc on the fitting curve. Then, the distance between 
data point Xfc and the fitting curve is ||P(£fc) — Xfc || . Here, 
tfc is called the location parameter of Xfc , P (tk ) is called the 
foot point corresponding to Xk. Denote T = {ti, ...,tk}, 
i.e. the collection of the location parameters of all the data 
points. The fitting problem is then formulated as: 
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where Ff airing is a fairing term which defines the fairness 
of a curve. Ff airing is commonly defined as follows [21]: 



P'(t)fdt + f3 / ||P"(t)|| 2 dt (2) 



fairing 



Since the objective function in Eqn. 1 is nonlinear, it is 
natural to apply iterative minimization methods to solve 
it. Most prevailing methods for solving this problem in 
CAGD are not standard optimization methods in the sense 
that they separately optimize location parameters T and 
control points V , making the problem much simpler to han- 
dle. However, these methods are time-consuming because 
they need to compute foot points on the fitting curve and 



to formulate and solve linear systems in every iteration. 
We observe that these time-consuming operations can be 
avoided by employing a L-BFGS optimization method that 
solves T and V simultaneously. We show that the result- 
ing algorithm is very efficient because in every iteration it 
does not need to perform foot point projection or solve a 
linear system of equations. 

The remainder of this paper is organized as follows. In 
section 2, we review some previous work. Section 3 intro- 
duces the standard L-BFGS optimization method. Section 
4 presents our new algorithm. Section 5 shows experimen- 
tal results and comparisons with existing methods. Then 
we conclude the paper in Section 6 with discussions of fu- 
ture work. 

2 Related Work 

Problem 1 can also be formulated as a nonlinear con- 
strained minimization problem with unknown variables V: 
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mn - 22 ||P(*fc) - Xfc|| 2 + \Ff airing, 



(3) 



where each tk is chosen such that P(tfc) is the foot point 
of Xfc and thus satisfies: 



(p(t fe )-x fe ) T p;(t fe ) = 0, k = l, 



(4) 



By representing {t k } as functions of V and putting them 
into the objective function in Eqn. 3, we obtain a nonlinear 
unconstrained minimization problem of variables V: 
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mn 2 ^ \\P(t k (V)) - X fe || 2 + XF fairing . 



(5) 



This is the viewpoint taken in [12] that reveals inherent 
relationship between some traditional methods and stan- 
dard optimization techniques. Most methods for solving 
problem 5 deal with control points and foot points sepa- 
rately [9] [12]. Each iteration of these methods consists of 
the following two steps: 

Step 1: Foot point projection: Fixing the control 
points of the current fitting curve, compute the location 
parameters T = {tk} for the data points {Xfc} such that 
{P(tfc)} are the foot points of {Xfc} on the current fitting 
curve. This step preserves the orthogonality constraint in 
Eqn. 4. 

Step 2: Control point update: In this step, T is fixed 
and a quadratic function e k in terms of the control points 
V is used to approximate the nonlinear squared distance 
from a data point Xfc to the fitting curve. Then the control 
points V are computed by minimizing the quadratic objec- 
tive function Q(V) = J2k e k(V) + \F s air ing ■ Since both e fe 
and Ffairing are quadratic functions of V, this step entails 
solving the linear equations VQ('P) = 0. 

Depending on different quadratic approximations cho- 
sen for efc, there are mainly three kinds of existing op- 
timization methods for curve fitting. The first one is 
the Point Distance Minimization method, or PDM. This 
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method is widely used because of its simplicity. Refer- 
ences on PDM include (but are not limited to) [14], [15], [8] 
and [19] on curve fitting as well as [7], [5] and [6] on surface 
fitting. The error term used in PDM is defined by 

e PD ,k = ||P(P;^)-X fc || 2 . (6) 

Geometrically, this function defines the distance between 
a data point and a point on the fitting curve at a partic- 
ular parameter t%- Considering the fact that t k (V) is set 
to a constant, this definition is a poor approximation of 
the nonlinear distance in Eqn. 5. As pointed out in [3], 
from the viewpoint of optimization, PDM is an alternat- 
ing method and exhibits linear convergence rate. We will 
see in our experiments that PDM is the slowest among all 
the methods we have tested. 

The second method is called the tangent distance min- 
imization method (TDM) [4] which uses the error term 

e T D,k= [(P(P;t fc )-X fc ) T -N fc ] 2 , ( 7 ) 

where Nfc is the unit normal vector at point P(t k ) on the 
curve. 

The term eTD,k defines the distance between a data 
point Xfe and the tangent line at P(t k ). Although this is a 
fair approximation to the true squared distance near a flat 
part of curve, it is not accurate near high curvature regions 
since no curvature information is considered. As a result, 
TDM does not show stable performance near high curva- 
ture regions [21]. In fact, it has been pointed out in [21] 
that TDM is essentially Gauss- Newton minimization with- 
out step-size control, and regularization should be used to 
improve the stability of TDM. 

Applying the Levenberg-Marquardt regularization to 
TDM leads to a method called TDMLM [21]. Suppose 
the linear system for control points updating in TDM is 
At dm ■ V = t>T d m , where At dm is a matrix and \)tdm 
is a vector. In TDMLM, the control points V are computed 
by solving 

(At DM + ^1) ■ V = t>TDM- 

Empirically, jj, is set as p — ir ^^°"^ , where tr(ATDM) is 
the trace of Atdm, n the number of control points, and I 
the identity matrix. 

The third method, called the Squared Distance Mini- 
mization method or SDM [21], uses a curvature-based error 
term, which is a variant of the second order approximation 
to the true squared distance introduced in [17] [16]. This 
error term, called the SD error term, is defined by 

[ ^ [(P(V;t k )-X k ) T -T k ] 2 + 
esD, k =< + [{P(P;t k )-X k f -N k ] 2 ,iid<0, (8) 
[ [{P(P;t k )-X k ) T -N k ] 2 ,ifO<d<p, 

where p is the curvature radius at P(t k ) and d is the pos- 
itive distance between Xfe and P(t k ). The SD error term 
contains some second order derivative information and is 
therefore a better approximation to the true squared dis- 
tance function than those used in TDM and PDM. From 
the viewpoint of optimization, SDM is quasi-Newton opti- 
mization method that employs a modified Hessian matrix 
of the original nonlinear distance function. This modifi- 
cation discards some complicated parts in the true Hes- 
sian matrix and keeps other parts with intuitive geomet- 
ric meanings [21]. The semi-definite positive property of 
the modified Hessian matrix is also guaranteed. It has 
been demonstrated in [21] that SDM exhibits better per- 
formance in terms of convergence rate and stability than 
PDM and TDM. 



Since the curve fitting problem is formulated as a non- 
linear least squares minimization problem in Eqn. 1, it is 
natural to study how to solve it using standard optimiza- 
tion methods. The authors of [12] apply the Gauss-Newton 
method to Eqn. 1 and derive new error terms using simpli- 
fied partial derivatives of the objective function in Eqn. 1. 
These methods are observed to have similar performances 
as SDM. 

All the above methods update control points V and 
location parameters {t k } in two interleaving steps. The 
main difference of our new method with these existing 
methods is that in every iteration we update V and {t k } 
simultaneously. In this sense the most closely related work 
is [20] which also optimizes control points and location pa- 
rameters simultaneously in every iteration. However, that 
method uses the Gauss-Newton optimization and therefore 
still needs to valuate and store the Jacobian matrices of the 
objective function, whose size depends on the number of 
data points and control points [20], as well as to solve a lin- 
ear system of equations. In contrast, our approach based 
on L-BFGS does not need to formulate and solve any linear 
equations and is therefore faster than the method in [20], 
as we are going to demonstrate in later experiments. 

Other optimization techniques have been explored for 
surface and curve fitting problems in literature. The au- 
thors of [22] proposed a method for NURBS curve and 
surface fitting which optimizes control points, parameters 
and knots by a conjugate gradient method. Genetic Algo- 
rithms and optimal control methods have also been tried 
in curve fitting [18] [2]. These methods are generally slow 
and have only been applied to simple examples. 

3 Limited Memory BFGS - L-BFGS 

Limited Memory BFGS, or L-BFGS, is a quasi-Newton 
method for solving unconstrained nonlinear minimization 
problems [13]. L-BFGS approximates the inverse Hessian 
matrix of the objective function by a sequence of gradient 
vectors from previous iterations. Suppose we want to solve 
an unconstrained optimization problem 

min/(a:), 

where f(x) is a nonlinear function to minimize and x a set 
of unknown variables. In the fc-th iteration of L-BFGS, the 
variables x k +i are updated by 

x k+1 = x k — a k H k \7f(x k ), 

where H k is an approximation to the inverse Hessian ma- 
trix of f(x) at x k . Here, —H k X7f(x k ) is a search direction, 
and Qffe a scalar variable controlling the step-size of search 
direction [13]. 

Define s k := x k+1 -x k , y k := \7f k+1 -\7f k , p k = 

V k = I — p k y k s k . L-BFGS uses the values of the objective 
function and its gradient in the (k— m)-th iteration through 
(k — l)-th iteration to compute H k [13]: 

H k = (Vf_! ■ ■ • V k T _ m )H° k (V k „ m ■ • • V h -i) 

+ Pk- m (V k _ 1 ■ ■ ■ V fc _ m+1 )Sfe_ m Sfe_ m (Vfe_ m + l • ■ • Vfe_l) 
+ p k -m + l(V k _ 1 ■ ■ ■ Vfe_ m + 2 )Sfc-m+l- 
•Sfe_ m + l(Vfe_ m + 2 ' ■ ' Vfe_l) 
+ ••• 

+Pfc-lSfe-lSfc-l, 

(9) 

where is a diagonal matrix defined by — 7*./, where 

T 

s k-\yk-l r-, ol 

7fc - vl^l [13] - 

In practice, we do not need to compute and store the 
matrix H k . Instead, we compute the search direction 
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— i/fcV/fc directly by a L-BFGS two-loop recursion algo- 
rithm (Algorithm 1) [13]: 



3. Run Algorithm 2: the L-BFGS algorithm until con- 
vergence. 



Algorithm 1 L-BFGS two-loop recursion 
q = V/(ar k ); 

for i — k — 1, k — 2, . . . , k — m do 

on = Pisjq; 
q = q- onyi; 
end for 

Z = H° k q; 

for i — k — m,k — m + 1, . . . k — 1 do 

z — z + Si(at — pi); 
end for 

Output z to be //fcV/fe. 



The L-BFGS optimization procedure is described in 
Algorithm 2 [13]. 



Algorithm 2 the L-BFGS algorithm 

Choose a starting point xo and a positive integer m; 

k = 0; 

repeat 

Choose fl£; 

Compute a descending direction pt by the two-loop 
recursion algorithm; 

Compute Xk+i = Xh + otkPk, otk is chosen to satisfy 
the Wolfe conditions; 
if fc > m then 

Discard s k - m and yk-m,; 
end if 

Compute the values of Sk, yt and store them; 
until convergence 



L-BFGS stops when the norm of the gradient of the ob- 
jective function is smaller than a specified tolerance value 
e, i.e. ||V/|| < e. 

In algorithm 2, once the descending direction pk is ob- 
tained, the variables Xk+i should be updated by Xk+i = 
Xk + ctkPk- Here is chosen to guarantee the decreas- 
ing of the value of the objective function. This is usually 
solved by a linesearch algorithm. Basically, a linesearch 
algorithm starts with otk = 1 and decreases the value of 
a.k by some strategy until the following Wolfe conditions 
are satisfied [13]: 

f f(x k + QfeP*,) < f(x k ) + c 1 a k Vf k r pk, 
1 V/(> fc + a k Pk) T Pk > c 2 V/Jp ft . 

Here ci and C2 are constants which satisfy < ci < C2 < 1. 
In our algorithm we use c\ = 10" 4 and c 2 = 0.9 through- 
out optimization. We have found from our experience that 
ctk = 1 is often good enough and the computational time 
of linesearch only takes a small partition of the total com- 
putational time. We will show these timing data in later 
sections. 

4 Curve fitting with L-BFGS 
4.1 Algorithm outline 

We employ L-BFGS directly to solve the nonlinear least 
squares minimization problem 1 in the following steps. 

1. Specify an initial curve P(t). 

2. Find the foot point P(ifc) on P(t) for every data point 
Xfc. This gives the initial position of location param- 
eters {tfc}. 



We will refer to this algorithm as the L-BFGS fitting 
method in the following sections. 

4.2 Selection of m 

The L-BFGS algorithm (Algorithm 2) uses m gradient vec- 
tors in a sequence of iterations to approximate the inverse 
Hessian matrix. A larger m can result in a more accu- 
rate approximation but at the same time, it will take more 
computational time. Therefore, it is important to select a 
proper value of m which balances between the competing 
objectives, i.e., accuracy and efficiency. In literature, the 
value of m is often chosen between 3 and 20 [13]. How- 
ever there are also papers reporting that a very large value 
of m is necessary for generating satisfactory results. One 
example is [11] in which m is set to 240. 

To find a proper value of m in our curve fitting prob- 
lem, we have tested many examples with m — 3, 5, 7, 20, 
50 and 120 to understand the behavior of our algorithm. 
Our conclusion is, for simple examples, using m — 3 and 
m = 5 may lead to a slightly faster fitting error descending 
speed. On the other hand, using m = 20, 50 and 120 will 
contribute to a faster gradient norm convergence speed. 
For more complicated examples (examples with more data 
points and heavier noise), the behaviors of error descending 
speed and gradient norm descending speed using different 
m tend to be indistinguishable. Considering all these fac- 
tors, we suggest using m = 20 in experiments. 

4.3 Foot point projection 

Computing foot points needs to be performed in every it- 
eration of traditional curve fitting methods. Foot point 
projection is itself an optimization problem which is inves- 
tigated in [9], [19], [1] and [10]: 

min||P(t)-X fe || 2 . 

When implementing the traditional methods we compute 
foot points in every iteration using the Gauss-Newton 
method outlined below. 

1. Initialization: In this step location parameters cor- 
responding to data points are roughly estimated. The es- 
timation will be used as initialization for further optimiza- 
tion in the next step. A straightforward method is: Sam- 
pling dense enough points on the fitting curve and finding 
the closest sample to each data point as the estimation of 
foot points. 

2. Iterative Update: Update the location parameter it- 
eratively [21]: in the i-th iteration, the location parameter 
of Xfe is updated by tk,i+i ~ tk,i + a8t, where 

(X fc -P(t M )).P'(t M ) 

I|p'(*m)II 2 

is the suggested update in the descending direction. The 
value of a is decided by a simple linesearch method to 
guarantee the decreasing of the orthogonal distance, i.e., 
||P(*M+i) _ x fc|l < l|P(*fc,0 — x fc|l- Tne optimization 
process stops when ||(X fe - P(tfc,i)) ' p '(*fe,i)ll < 10~ 10 . 

The Initialization step is generally time consuming. In 
practice, after several iterations in the beginning of curve 
fitting, the shape of the fitting curve will not change a 
lot in later optimization. In this case, we can use foot 
points computed in the i-th iteration (i.e.,{tfc,i}) as initial- 
ization foot points for the (i-fl)-th iteration. To determine 
whether it is safe to directly use the foot points from the 
last iteration as initialization, we use a criterion based on 
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variation of fitting error: We suggest to re-initialize foot 
points in the (i + f)-th iteration when — t ^r~~ > 0.2 
where Ei and Ei+i are the fitting error in the i-th iteration 
and (i + l)-th iteration respectively, as defined in Eqn. 10. 
In our experiments, with this criterion, for all traditional 
methods, the number of foot point initializations needed is 
from 1 to 4 in all our tested examples. 

4.4 Foot point correction and restart of L-BFGS 

The initialization of the L-BFGS fitting method also needs 
foot point projection to determine the initial set T of lo- 
cation parameters. Then, in the subsequent iterations the 
L-BFGS fitting method in general does not need to perform 
foot point projection any more but rather optimizes loca- 
tion parameters {tk} and control points V simultaneously. 
In rare case, especially when the initial fitting curve speci- 
fied is not good enough, it is possible that P(ife) is far from 
the closet point on the curve to Xfc, even if P(tfc) — Xfc is 
orthogonal to P'(tfc). An example is shown in Figure 1(a): 
it is part of a fitting curve on convergence, but there are 
points P(tfc) which are not the foot points of Xfc, because 
the L-BFGS fitting method gets stuck in a poor local min- 
imum. In this case the following remedy can be used. 

To rectify the incorrect foot point projections, we just 
perform foot point computation after the termination of 
current L-BFGS algorithm, and start a new L-BFGS algo- 
rithm taking initial control points from the previous run of 
L-BFGS and initial location parameters from the output of 
foot point computation. Figure 1(b) shows the result fit- 
ting curve. To detect a local minimum automatically, we 
measure the fitting error E after the L-BFGS algorithm 
and the fitting error E + after the additional operation of 
foot point computation. If the error changing \\E — E+\\ is 
bigger than a tolerance (we use 10 -6 in our experiments), 
we conclude that the foot points of some data points are 
corrected and an additional run of the L-BFGS algorithm 
is needed. 





(a) Without foot 
point correction. 



(b) With foot 
point correction. 



Figure 1: foot point correction. 

5 Results and Discussions 

In this section, we first present some experiments com- 
paring the L-BFGS fitting method with existing methods, 
then we give explanation on the fast speed of the L-BFGS 
fitting method. 

5.1 Experiments 

The fitting error of the i-th iteration is measured by: 



E, 



N 1 

£-||P(t fc )-Xfc 



(10) 



The parameter domains of B-spline curves in these ex- 
amples are set to [0, 1]. All data points are scaled into a 
unit box: [0, 1] x [0, 1]. 

Due to different complexities and the set up of initial 
curves in the examples in our experiments, we use differ- 
ent coefficients of fairing terms a and f3 in different exam- 
ples to obtain satisfactory fitting curves. In each example, 
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(b) Fitting curve. 
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(g) Time cost for an iteration (in seconds). 

Figure 2: The target shape is a set of 100 points on a circle. 
A B-spline curve with 6 control points is used to fit it. No 
fairing term is used in this example. 

the same values of faring term coefficients are used for all 
tested methods. The values of coefficients are noted in the 
captions of Figures. 

Comparison with traditional methods. Three 
data sets are given in Figure 2, 3 and 4 for comparisons 
with three traditional methods: PDM, TDMLM and SDM. 
For each data set, we show data points, the initial fit- 
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(b) Fitting curve. 



(a) Initialization. 



(b) Fitting curve. 



0.03 
0.02 




TDMLM 

SDM 

PDM 

L-BFGS 



100 200 300 400 500 600 700 

iterations 

(c) Error vs iteration. 



0.03 
0.02 



TDMLM 
SDM 
PDM 
L-BFGS 



5 10 15 20 25 30 35 40 45 50 

iterations 

(c) Error vs iteration. 



0.03 
0.02 



1e + 
1e-1 
1e-2 
1e-3 
1e-4 
1e-5 
1e-6 
1e-7 
1o-S 




0.0 0.2 0.4 0.6 0.8 

time(s) 
(d) Error vs time. 



1.0 1.2 1.4 



TDMLM 
SDM 
PDM 
L-BFGS 




0.0 0.2 0.4 0.6 0.8 

time(s) 

(e) Gradient norm vs time. 
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(g) Time cost for an iteration (in seconds). 

Figure 3: An example with sharp features. The fitting 
curve has 24 control points and the data set contains 90 
points. The coefficients of fairing term are set to a = 
and 8 = 5 ■ 10~ 4 . 
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(g) Time cost for an iteration (in seconds). 

Figure 4: An example with noisy data points. The fitting 
curve has 8 control points and the data set has 150 points. 
No fairing term is used in this example. 



ting curve and the final fitting curve of the L-BFGS fitting 
method. Three charts are also shown for each data set. 
The first two charts show the fitting error versus the it- 
eration number and computational time respectively. The 
third chart shows the decreasing of gradient norm versus 
computational time. 



We observe that in the first several iterations, the fit- 
ting error of the L-BFGS fitting method does not decrease 
as fast as SDM and TDMLM in terms of number of it- 
erations. That is because that in the L-BFGS algorithm 
(Algorithm 2), the approximation of inverse Hessian ma- 
trix needs to be accumulated by using information from 
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a sequence of m iterations. Therefore, at the first several 
iterations, the approximant matrix is not accurate enough 
and this slows down the performance of the L-BFGS fit- 
ting method. However, an iteration of the L-BFGS fitting 
method is much faster than PDM, TDMLM and SDM. As 
a result, the L-BFGS fitting method converges much faster 
than the other three methods in terms of computational 
time, as shown in Figure 2(d), 3(d) and 4(d). 

Convergence. The convergence behaviors of the four 
methods can be observed from the third chart in the above 
three examples, showing the decrease of gradient norm 
against computational time. The termination criterion for 
all these examples is || V/||oo < 1CF 8 , where / is the objec- 
tive function. From Figure 2(e), 3(e) and 4(e), we observe 
that the L-BFGS fitting method is the only method that 
always meets this criterion, i.e. the gradient norm of its ob- 
jective function reaches the threshold of 10 -8 . This is not 
surprising since the L-BFGS fitting method implements a 
well-studied optimization method (the L-BFGS algorithm, 
Algorithm 2) that has demonstrated superior convergence 
behavior close to the superlinear rate possessed by the 
BFGS method [13]. 




(a) Initializa- (b) Result of L- (c) Result of (d) Result of 
tion. BFGS. SDM. Speer's method. 





L-BFGS SDM 


Speer's 


time(s) 
error 


1.0202 9.2754 
0.01796 0.01805 


87.5810 
0.01947 



(e) Timing and fitting errors. 



Figure 5: Comparisons with L-BFGS, SDM and Speer's 
method. The coefficients of fairing term are a = and 
/3 = 10" 3 . 

Comparison with the method of Speer et al. 

In [20], Speer et al proposed to use the Gauss-Newton 
method to solve the least squares problem 1. This method 
also optimizes T and V simultaneously. However, in every 
iteration it still needs to formulate and solve linear sys- 
tem which includes both location parameters and control 
points as variables. Therefore, this method is inefficient 
for large data sets. Figure 5 shows an example with noise 
containing 2500 data points. We use this example to com- 
pare the L-BFGS fitting method with SDM and Speer's 
method. From Table 5 we can see that the L-BFGS fitting 
method is capable of producing a satisfactory curve about 
9 times faster than SDM; SDM in turn is about 8 times 
faster than Speer's method. 

More examples. We present more examples in Fig- 
ure 8 and 9. 

5.2 Analysis and discussions 

It is difficult to provide a theoretical proof on the superior 
efficiency of the L-BFGS fitting method over existing meth- 
ods. As an alternative, in this section we shall conduct 
an empirical study on the efficiency of PDM, TDMLM, 
SDM and the L-BFGS fitting method, in order to gain a 
better understanding of their relative performances. The 
traditional methods (PDM, TDMLM, SDM) that update 
control points and location parameters separately mainly 
include the following tasks: linear system formulation and 
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Matrix 


foot point 






filling 


solving 


projection 


Example 2 


PDM 


23.8% 


15.5% 


60.6% 


(ctrl pts: 6 


TDMLM 


22.6% 


35.3% 


41.9% 


data pts: 100) 


SDM 


31.5% 


41.4% 


26.8% 


Example 3 


PDM 


10.7% 


55.6% 


33.8% 


(ctrl pts: 24 


TDMLM 


7.1% 


67.3% 


26.0% 


data pts: 90) 


SDM 


8.8% 


70.0% 


21.1% 


Example 4 


PDM 


25.5% 


26.6% 


48.3% 


(ctrl pts: 8 


TDMLM 


25.5% 


29.6% 


45.2% 


data pts: 150) 


SDM 


33.9% 


26.6% 


39.8% 


Example 8 


PDM 


21.2% 


34.2% 


45.0% 


(ctrl pts: 30 


TDMLM 


8.1% 


78.6% 


13.0% 


data pts: 600) 


SDM 


23.1% 


50.2% 


27.2% 


Example 9 


PDM 


9.9% 


15.1% 


73.3% 


(ctrl pts: 66 


TDMLM 


15.6% 


31.8% 


50.0% 


data pts: 2000) 


SDM 


29.0% 


26.9% 


41.9% 



Table 1: Computational time for different parts of the PDM, 
TDMLM and SDM methods. 





Computing des- 
cending direction 


Linesearch 


Example 2 


95.8% 


4.2% 


Example 3 


97.8% 


2.1% 


Example 4 


96.8% 


3.2% 


Example 8 


91.4% 


8.6% 


Example 9 


89.7% 


10.3% 



Table 2: Computational time for different parts of the L-BFGS 
algorithm. m=20. 

solving for control points and foot points computation for 
location parameters. The timing data for different parts of 
the methods for the examples in this paper are presented 
in Table 1. The L-BFGS fitting method consists of two 
parts: the two-loop algorithm for computing a descending 
direction and a line-search algorithm for deciding step-size. 
Timings for these two parts on the same examples as in Ta- 
ble 1 are listed in Table 2. 

We have the following observations on these timing 
data. 

• Although the number of control points is generally 
much fewer than the number of data points, tradi- 
tional methods still consume more than half of the 
total time on updating control points, because of the 
need to fill the matrix and solving the linear sys- 
tem in every iteration. The L-BFGS algorithm (Al- 
gorithm 2) is a Newton-type optimization method 
that uses an approximated inverse Hessian matrix of 
the objective function. However, instead of solving a 
large linear system to compute the descend direction 
as PDM, TDMLM and SDM, the L-BFGS algorithm 
uses a two-loop algorithm which uses only vector mul- 
tiplications and is therefore much faster. 

• Foot point computation is very time consuming. If 
we re-compute the initialization of foot points in ev- 
ery iteration, the overall time for foot point computa- 
tion would be more than 90% of the total time of the 
algorithm, as observed in [21]. In our implementa- 
tion of the traditional methods used for comparison 
in this paper, we use as much as possible the foot 
points in the previous iteration as initialization for 
the current iteration, thus having saved a lot of time 
for traditional methods. Even so, the L-BFGS fit- 
ting methods still outperforms these traditional meth- 
ods, since there is generally no need to perform foot 
point projection in the L-BFGS fitting method. In 
rare cases, foot point computation is needed for the 



6 



L-BFGS method to jump out of a poor local mini- 
mum, as we have explained in section 4. This is an 
issue mostly due to the quality of initialization, rather 
than the inherent demand of the algorithm. 

• The L-BFGS fitting method performs optimization 
in a much higher dimensional space than those of tra- 
ditional methods since generally the number of data 
points is much larger than that of the control points. 
Therefore, the terrain of the functional is supposed 
to be much more complicated and the optimization 
is more difficult. The linesearch algorithm is there- 
fore necessary for stable convergence of the L-BFGS 
fitting algorithm. Table 2 shows that the computa- 
tional time by the linesearch algorithm usually takes 
a small part of the total time (less than 10% in most 
cases). 

We now study how the computational time depends on 
the number of control points and and the number of data 
points. 

Timing vs # of data points. In Figure 6, we show 
computational time with increased number of data points 
for various methods. The number of data points in these 
point sets is 100, 200, 500, 1000 and 3000 respectively. 
The fitting curve has 8 control points. Figure 6 shows that 
the computational time for each iteration of all 4 methods 
depends almost linearly on the number of data points. This 
can be explained as follows. It is not difficult to see that in 
the PDM, TDMLM and SDM, the time for matrix building 
and foot point projection is linear in the number of data 
points. The time for solving linear system is constant since 
the number of control points is fixed. Consequently, the 
total time for these three methods increase linearly as the 
number of data points increases. For the L-BFGS fitting 
method, computational time is linear in the number of 
variables (2 x the number of control points + the number 
of data points), therefore the computational time of the L- 
BFGS fitting method also increases linearly as the number 
of data points increases. 

Timing vs # of control points. The relationship 
of computational time and the number of control points of 
the fitting B-spline curve can be observed in Figure 7. We 
insert new control points by knot insertion in each knot in- 
terval and get 4 B-spline curves with the number of control 
points: 10, 20, 40 and 80 respectively. The number of tar- 
get data points is 200. We see that the per-iteration time 
for the tested traditional methods increases faster than the 
L-BFGS fitting method when the number B-spline control 
points increases. That is because in the PDM, TDMLM 
and SDM, the size of linear system is quadratic to the 
number of control points, but the computational time of 
the L-BFGS algorithm (i.e. the two-loop algorithm and 
the linesearch) depends on the number of control points 
linearly. 

These experiments show that the L-BFGS fitting 
method is more suitable for large scale curve fitting prob- 
lems, especially when the target shape is complicated and 
a large number of control points are involved. 




(a) 100 data points, (b) 200 data points, (c) 500 data points. 





(d) 1000 data points. 



(c) 3000 data points. 
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(f) Per-iteration time as the number of data points increases. 

Figure 6: Increasing data points: An example with 8 con- 
trol points. The coefficients of fairing term are a = 5- 10 -4 
and = 




(a) 10 Con- 
trol points 



(b) 20 Con- 
trol points 



(c) 40 Con- 
trol points 



(d) 80 Con- 
trol points 





TDMLM 






SDM 




PDM 




L-BFGS 




20 40 60 80 

number of control points 

(e) Per-iteration time as the number of control points in- 
creases. 

Figure 7: With increased number of control points. There 
are 200 data points. No fairing term is used. 



6 Conclusion and Future Work 

In this paper, we propose a new curve fitting method based 
on the L-BFGS optimization technique. The unique fea- 
tures of this algorithm are that it does not need to perform 
the time-consuming foot point projection in every iteration 
as in traditional approaches and that it does not need to 
formulate and solve a linear system of equations in every 
iteration; instead, it uses only efficient vector multiplica- 
tions. As a result, this new method is much faster than 



other traditional methods, as demonstrated by a number 
of experimental results presented. In the future we will 
extend this method to solving the B-spline surface fitting 
problem, for which we expect even more significant im- 
provements over the existing methods because of the large 
number of data points as well as the large number of con- 
trol points involved in surface fitting. 
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(a) Initialization. 



(b) Fitting curve. 




TDMLM 

SDM 
PDM 
L-BFGS 



I 0.1 0.2 0.3 0.4 

time(s) 

(c) Fitting error vs time. 



L-BFGS PDM 



8.2- 10" 



5.2 



TDMLM 
0.21 



SDM 
0.23 



(d) Time to attain minimal error (in seconds). 

Figure 8: A Chinese character with 30 control points and 
600 data points which means "mountain". The coefficients 
of fairing term are a = 5 ■ 10~ 4 and ,3 = 0. 




(a) Initialization. 



0.03 
0.02 



(b) Fitting curve. 




TDMLM 
SDM 
PDM 
L-BFGS 



I 0.5 1.0 1.5 

time(s) 

(c) Fitting error vs time. 



L-BFGS PDM TDMLM SDM 
0.28 41 1.4 1.8 

(d) Time to attain minimal error (in seconds). 



Figure 9: Flame with 66 control points and 2000 data 
points. The coefficients of fairing term are a = 10~ 3 and 
13 = 10" 2 . 
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